#####################################################################
# Figure 1
#####################################################################

#####################################################################
# Preliminaries
#####################################################################

rm(list = ls())

library(scinference)
library(limSolve)

set.seed(12345)

###################################################################
# Functions
###################################################################

source("code/auxiliary scripts/common_functions.R")

sim <- function(T0,T1,J,K,rho.u,mu){
  
  # Preliminaries
  alpha.sig <- 0.1
  
  # DGP
  T01       <- T0+T1
  
  w <- rep(0,J)
  w[1:3] <- 1/3
  
  Y0        <- matrix(rnorm(T01*J),T01,J)
  Y0[,1:3]  <- Y0[,1:3] + 2 
  Y1        <- mu + Y0 %*% w + generate.AR.series(T01,rho.u)
  
  # SC without debiasing
  Y1.pre  <- Y1[1:T0]
  Y0.pre  <- Y0[1:T0,]
  Y1.post <- Y1[(T0+1):T01]
  Y0.post <- Y0[(T0+1):T01,]
  
  w.hat     <- scinference:::sc(Y1.pre,Y0.pre,lsei_type=1)$w.hat
  tau.nodb  <- mean(Y1.post -Y0.post %*% w.hat)
  
  # SC with debiasing
  tau.db <- scinference(Y1,Y0,T1=T1,T0=T0,inference_method="ttest",estimation_method="sc",K=K,lsei_type=1)$att
  
  return(c(tau.nodb,tau.db))
  
}


###################################################################
# Simulations
###################################################################

# Setup (based on empirical application)
nreps <- 50000
K     <- 3
T0    <- 30 
T1    <- 16 
J     <- 14 
rho.u <- 0.31

# Simulation loops
result_correct <- matrix(NA,nreps,2)
for (r in 1:nreps){
  result_correct[r,] <- sim(T0=T0,T1=T1,J,K,rho.u,0)
}

result_misspec <- matrix(NA,nreps,2)
for (r in 1:nreps){
  result_misspec[r,] <- sim(T0=T0,T1=T1,J,K,rho.u,2)
}

###################################################################
# Graphics
###################################################################

original_correct <- result_correct[,1]
debiased_correct <- result_correct[,2] 

original_misspec <- result_misspec[,1]
debiased_misspec <- result_misspec[,2] 

xl <- -2
xu <- 2

yl <- 0
yu <- 1.2

graphics.off()
pdf("graphics/hist_sc.pdf",pointsize=16,width=8.0,height=6.0)
h <- hist(original_correct, breaks=c(seq(xl,xu,length=40)), col="gray", xlab="", main="Correct specification: SC estimator w/o debiasing",freq=F,xlim=c(xl,xu),ylim=c(yl,yu))
xfit <- seq(min(xl),max(xu),length=1000) 
yfit <- dnorm(xfit,mean=0,sd=sd(original_correct)) 
lines(xfit, yfit, col="black", lwd=5)
text(x=1.5,y=1.1, label="bias =",adj=1)
text(x=1.75,y=1.1, label=sprintf("%.3f",mean(original_correct)))
text(x=1.5,y=1, label="sd =",adj=1)
text(x=1.75,y=1, label=sprintf("%.3f",sd(original_correct)))
legend("topleft",legend=c("Simulation","Normal curve"), col=c("gray","black"),fill=NA,border=NA, lty=c(1,1),lwd=c(15,5),merge=T,bty="n")
dev.off()

graphics.off()
pdf("graphics/hist_debiased.pdf",pointsize=16,width=8.0,height=6.0)
h <- hist(debiased_correct, breaks=c(seq(xl,xu,length=40)), col="gray", xlab="", main="Correct specification: Debiased SC estimator",freq=F,xlim=c(xl,xu),ylim=c(yl,yu)) 
xfit <- seq(min(xl),max(xu),length=1000) 
yfit <- dnorm(xfit,mean=0,sd=sd(debiased_correct)) 
lines(xfit, yfit, col="black", lwd=5)
text(x=1.5,y=1.15, label="bias =",adj=1)
text(x=1.75,y=1.15, label=sprintf("%.3f",mean(debiased_correct)))
text(x=1.5,y=1.05, label="sd =",adj=1)
text(x=1.75,y=1.05, label=sprintf("%.3f",sd(debiased_correct)))
legend("topleft",legend=c("Simulation","Normal curve"), col=c("gray","black"),fill=NA,border=NA, lty=c(1,1),lwd=c(15,5),merge=T,bty="n")
dev.off()

xl <- -2.5
xu <- 4.5

yl <- 0
yu <- 1.2

graphics.off()
pdf("graphics/hist_sc_misspec.pdf",pointsize=16,width=8.0,height=6.0)
h <- hist(original_misspec, breaks=c(seq(xl,xu,length=40)), col="gray", xlab="", main="Misspecification: SC estimator w/o debiasing",freq=F,xlim=c(xl,xu),ylim=c(yl,yu))
xfit <- seq(min(xl),max(xu),length=1000) 
yfit <- dnorm(xfit,mean=0,sd=sd(original_misspec)) 
lines(xfit, yfit, col="black", lwd=5)
text(x=3.6,y=1.15, label="bias =",adj=1)
text(x=4,y=1.15, label=sprintf("%.3f",mean(original_misspec)))
text(x=3.6,y=1.05, label="sd =",adj=1)
text(x=4,y=1.05, label=sprintf("%.3f",sd(original_misspec)))
legend("topleft",legend=c("Simulation","Normal curve"), col=c("gray","black"),fill=NA,border=NA, lty=c(1,1),lwd=c(10,3),merge=T,bty="n")
dev.off()

graphics.off()
pdf("graphics/hist_debiased_misspec.pdf",pointsize=16,width=8.0,height=6.0)
h <- hist(debiased_misspec, breaks=c(seq(xl,xu,length=40)), col="gray", xlab="", main="Misspecification: Debiased SC estimator",freq=F,xlim=c(xl,xu),ylim=c(yl,yu)) 
xfit <- seq(min(xl),max(xu),length=1000) 
yfit <- dnorm(xfit,mean=0,sd=sd(debiased_misspec)) 
lines(xfit, yfit, col="black", lwd=5)
text(x=3.6,y=1.15, label="bias =",adj=1)
text(x=4,y=1.15, label=sprintf("%.3f",mean(debiased_misspec)))
text(x=3.6,y=1.05, label="sd =",adj=1)
text(x=4,y=1.05, label=sprintf("%.3f",sd(debiased_misspec)))
legend("topleft",legend=c("Simulation","Normal curve"), col=c("gray","black"),fill=NA,border=NA, lty=c(1,1),lwd=c(10,3),merge=T,bty="n")
dev.off()

